SUBROUTINE DG_RADFLX
!-----------------------------------------------------------------------
!
! Definition von DG:  Strahlungsflussgleichung
!

      use primvar,  only : X, XZ, XA, DG, MR, MD, ME, MUr, MJ, MH, &
                           DG_i, DG_im1, DG_ip1, DG_im2, DG_ip2
      use physco,   only : z1, z2, z3, z4, z12, z14, clight, pi, sqrt2pi
      use config,   only : np
      use global,   only : tst, zz
      use geomvar,  only : V_vol, V_flux, S_vol
      use eddvar,   only : f_edd, mu_r
      use zvar,     only : rho0, H_p, dH_pdE, dH_pdR
      use advecvar, only : H_adv, dH_adv_dim1, dH_adv_di, dH_adv_dip1, dH_adv_dip2
      use matvar,   only : OPAros, dOPArosdE, dOPArosdD
      use RBvar,    only : L_star


      implicit none

      integer :: i


!-----------------------------------------------------------------------
!    Randbedingungen
!-----------------------------------------------------------------------

   ! innere Pseudozellen: i=1 & i=2
      DG(MH,MH,DG_i,1)      =  z1
      DG(MH,MH,DG_ip1,1)    = -z1
      DG(MH,MH,DG_i,2)      =  z1
      DG(MH,MH,DG_ip1,2)    = -z1

   ! aeussere Pseudozellen: i=np & i=np-1
      DG(MH,MH,DG_i,np)     =  z1
      DG(MH,MH,DG_im1,np)   = -z1
      DG(MH,MH,DG_i,np-1)   =  z1
      DG(MH,MH,DG_im1,np-1) = -z1

   ! innere Randbedingung: i=3
      DG(MH,MH,DG_i,3)      =  z1
      !~ DG(MH,MR,DG_i,i)   =  z2 * L_star / (z4*pi)**2 / X(MR,i)**3 ! Einstrahlung

   ! aeussere Randbedingung: i=np-2
      DG(MH,MH,DG_i,np-2)   =  z1

!~ !-----------------------------------------------------------------------
!~ !    Mathematica Output: Abstrahlungsrandbedingung
!~ !-----------------------------------------------------------------------
!~ ! .... DG_radflx_RB.nb   mit:  G = Hi - z12 (XSvol Ji + XSvolim1 Jim1) / XVvol mueq
!~ !
!~ DG(MH,MR,DG_im1,i)=-((mu_r*(X(MJ,i)*S_vol(i) + X(MJ,i-1)*S_vol(i-1))) / (pi*(X(MR,i-1) + X(MR,i+1))*(X(MR,i+1) &
!~ - X(MR,i-1))**2)) + (mu_r*(X(MJ,i)*S_vol(i) + X(MJ,i-1)*S_vol(i-1))) / (pi*(X(MR,i-1) + X(MR,i+1))**2*(X(MR,i+1) - X(MR,i-1))) &
!~ - (mu_r*(-(X(MJ,i-1)*pi*(X(MR,i) + X(MR,i-1))) + X(MJ,i-1)*pi*(X(MR,i) - X(MR,i-1)))) / (pi*(X(MR,i-1) + X(MR,i+1))*(X(MR,i+1) &
!~ - X(MR,i-1)))
!~ DG(MH,MR,DG_i,i)=-((mu_r*(X(MJ,i-1)*pi*(X(MR,i) + X(MR,i-1)) - X(MJ,i)*pi*(X(MR,i) + X(MR,i+1)) + &
!~ X(MJ,i)*pi*(X(MR,i+1) - X(MR,i)) + X(MJ,i-1)*pi*(X(MR,i) - X(MR,i-1)))) / (pi*(X(MR,i-1) + X(MR,i+1))*(X(MR,i+1) - X(MR,i-1))))
!~ DG(MH,MR,DG_ip1,i)=(mu_r*(X(MJ,i)*S_vol(i) + X(MJ,i-1)*S_vol(i-1))) / (pi*(X(MR,i-1) + X(MR,i+1))*(X(MR,i+1) - &
!~ X(MR,i-1))**2) + (mu_r*(X(MJ,i)*S_vol(i) + X(MJ,i-1)*S_vol(i-1))) / (pi*(X(MR,i-1) + X(MR,i+1))**2*(X(MR,i+1) - X(MR,i-1))) - &
!~ (mu_r*(X(MJ,i)*pi*(X(MR,i) + X(MR,i+1)) + X(MJ,i)*pi*(X(MR,i+1) - X(MR,i)))) / (pi*(X(MR,i-1) + X(MR,i+1))*(X(MR,i+1) - &
!~ X(MR,i-1)))
!~ DG(MH,MJ,DG_im1,i)=-((mu_r*(X(MR,i) + X(MR,i-1))*(X(MR,i) - X(MR,i-1))) / ((X(MR,i-1) + X(MR,i+1))*(X(MR,i+1) &
!~ - X(MR,i-1))))
!~ DG(MH,MJ,DG_i,i)=-((mu_r*(X(MR,i) + X(MR,i+1))*(X(MR,i+1) - X(MR,i))) / ((X(MR,i-1) + X(MR,i+1))*(X(MR,i+1) - &
!~ X(MR,i-1))))
!~ DG(MH,MH,DG_i,i)=z1


!-----------------------------------------------------------------------
!    Restlicher Bereich
!-----------------------------------------------------------------------

      do i=4,np-3

DG(MH,MR,DG_im1,i)=-(X(MH,i)*pi*(X(MR,i-1) + X(MR,i+1))*z12 ) + X(MH,i)*pi*(X(MR,i+1) - X(MR,i-1))*z12  + &
(clight*tst*z14 *(-(X(MJ,i-1)*pi*(X(MR,i) + X(MR,i-1))*(-z1 **2 + f_edd(i-1)*z3 )) + X(MJ,i-1)*pi*(X(MR,i) - X(MR,i-1))*(-z1 **2 &
+ f_edd(i-1)*z3 ))) / X(MR,i) - (-(pi*X(MR,i-1)*z12 *z2 ) + pi*tst*XZ(MUr,i-1)*zz )*H_adv(i-1) + clight*X(MH,i)*tst*z12 &
*(-(pi*rho0(i-1)*(X(MR,i) + X(MR,i-1))*OPAros(i-1)) + pi*rho0(i-1)*(X(MR,i) - X(MR,i-1))*OPAros(i-1) - &
(X(MD,i-1)*S_vol(i-1)*OPAros(i-1)*dH_pdR(i-1)) / (sqrt2pi*H_p(i-1)**2) - (X(MD,i-1)**2*S_vol(i-1)*dH_pdR(i-1)*dOPArosdD(i-1)) / &
(sqrt2pi**2*H_p(i-1)**3))
DG(MH,MR,DG_i,i)=X(MH,i)*pi*tst*(X(MUr,i+1) - X(MUr,i-1)) + clight*pi*tst*(f_edd(i)*X(MJ,i) - &
f_edd(i-1)*X(MJ,i-1))*z2  - (clight*tst*z14 *(X(MJ,i)*S_vol(i)*(-z1 **2 + f_edd(i)*z3 ) + X(MJ,i-1)*S_vol(i-1)*(-z1 **2 + &
f_edd(i-1)*z3 ))) / X(MR,i)**2 + (clight*tst*z14 *(-(X(MJ,i)*pi*(X(MR,i) + X(MR,i+1))*(-z1 **2 + f_edd(i)*z3 )) + &
X(MJ,i)*pi*(X(MR,i+1) - X(MR,i))*(-z1 **2 + f_edd(i)*z3 ) + X(MJ,i-1)*pi*(X(MR,i) + X(MR,i-1))*(-z1 **2 + f_edd(i-1)*z3 ) + &
X(MJ,i-1)*pi*(X(MR,i) - X(MR,i-1))*(-z1 **2 + f_edd(i-1)*z3 ))) / X(MR,i) + (-(pi*X(MR,i)*z12 *z2 ) + pi*tst*XZ(MUr,i)*zz &
)*H_adv(i) - (-(pi*X(MR,i)*z12 *z2 ) + pi*tst*XZ(MUr,i)*zz )*H_adv(i-1) + clight*X(MH,i)*tst*z12 *(-(pi*rho0(i)*(X(MR,i) + &
X(MR,i+1))*OPAros(i)) + pi*rho0(i)*(X(MR,i+1) - X(MR,i))*OPAros(i) + pi*rho0(i-1)*(X(MR,i) + X(MR,i-1))*OPAros(i-1) + &
pi*rho0(i-1)*(X(MR,i) - X(MR,i-1))*OPAros(i-1) - (X(MD,i)*S_vol(i)*OPAros(i)*dH_pdR(i)) / (sqrt2pi*H_p(i)**2) - &
(X(MD,i)**2*S_vol(i)*dH_pdR(i)*dOPArosdD(i)) / (sqrt2pi**2*H_p(i)**3))
DG(MH,MR,DG_ip1,i)=X(MH,i)*pi*(X(MR,i-1) + X(MR,i+1))*z12  + X(MH,i)*pi*(X(MR,i+1) - X(MR,i-1))*z12  + &
(clight*tst*z14 *(X(MJ,i)*pi*(X(MR,i) + X(MR,i+1))*(-z1 **2 + f_edd(i)*z3 ) + X(MJ,i)*pi*(X(MR,i+1) - X(MR,i))*(-z1 **2 + &
f_edd(i)*z3 ))) / X(MR,i) + (-(pi*X(MR,i+1)*z12 *z2 ) + pi*tst*XZ(MUr,i+1)*zz )*H_adv(i) + clight*X(MH,i)*tst*z12 &
*(pi*rho0(i)*(X(MR,i) + X(MR,i+1))*OPAros(i) + pi*rho0(i)*(X(MR,i+1) - X(MR,i))*OPAros(i))
DG(MH,MD,DG_im1,i)=clight*X(MH,i)*tst*z12 *((S_vol(i-1)*OPAros(i-1)) / (sqrt2pi*H_p(i-1)) + &
(X(MD,i-1)*S_vol(i-1)*dOPArosdD(i-1)) / (sqrt2pi**2*H_p(i-1)**2))
DG(MH,MD,DG_i,i)=clight*X(MH,i)*tst*z12 *((S_vol(i)*OPAros(i)) / (sqrt2pi*H_p(i)) + &
(X(MD,i)*S_vol(i)*dOPArosdD(i)) / (sqrt2pi**2*H_p(i)**2))
DG(MH,ME,DG_im1,i)=clight*X(MH,i)*tst*z12 *(-((X(MD,i-1)*S_vol(i-1)*OPAros(i-1)*dH_pdE(i-1)) / &
(sqrt2pi*H_p(i-1)**2)) + rho0(i-1)*S_vol(i-1)*(-((X(MD,i-1)*dOPArosdD(i-1)*dH_pdE(i-1)) / (sqrt2pi*H_p(i-1)**2)) + &
dOPArosdE(i-1)))
DG(MH,ME,DG_i,i)=clight*X(MH,i)*tst*z12 *(-((X(MD,i)*S_vol(i)*OPAros(i)*dH_pdE(i)) / (sqrt2pi*H_p(i)**2)) + &
rho0(i)*S_vol(i)*(-((X(MD,i)*dOPArosdD(i)*dH_pdE(i)) / (sqrt2pi*H_p(i)**2)) + dOPArosdE(i)))
DG(MH,MUr,DG_im1,i)=-(X(MH,i)*pi*X(MR,i)*tst) - pi*XZ(MR,i-1)*tst*zz *H_adv(i-1)
DG(MH,MUr,DG_i,i)=pi*XZ(MR,i)*tst*zz *H_adv(i) - pi*XZ(MR,i)*tst*zz *H_adv(i-1)
DG(MH,MUr,DG_ip1,i)=X(MH,i)*pi*X(MR,i)*tst + pi*XZ(MR,i+1)*tst*zz *H_adv(i)
DG(MH,MJ,DG_im1,i)=-(clight*f_edd(i-1)*pi*X(MR,i)*tst*z2 ) + (clight*S_vol(i-1)*tst*z14 *(-z1 **2 + &
f_edd(i-1)*z3 )) / X(MR,i)
DG(MH,MJ,DG_i,i)=clight*f_edd(i)*pi*X(MR,i)*tst*z2  + (clight*S_vol(i)*tst*z14 *(-z1 **2 + f_edd(i)*z3 )) / &
X(MR,i)
DG(MH,MH,DG_im2,i)=-(V_flux(i-1)*zz *dH_adv_dim1(i-1))
DG(MH,MH,DG_im1,i)=-(V_flux(i-1)*zz *dH_adv_di(i-1)) + V_flux(i)*zz *dH_adv_dim1(i)
DG(MH,MH,DG_i,i)=V_vol(i) + pi*X(MR,i)*tst*(X(MUr,i+1) - X(MUr,i-1)) + clight*tst*z12 &
*(rho0(i)*S_vol(i)*OPAros(i) + rho0(i-1)*S_vol(i-1)*OPAros(i-1)) - V_flux(i-1)*zz *dH_adv_dip1(i-1) + V_flux(i)*zz *dH_adv_di(i)
DG(MH,MH,DG_ip1,i)=-(V_flux(i-1)*zz *dH_adv_dip2(i-1)) + V_flux(i)*zz *dH_adv_dip1(i)
DG(MH,MH,DG_ip2,i)=V_flux(i)*zz *dH_adv_dip2(i)

      end do


END SUBROUTINE DG_RADFLX




